Uniform Shock Waves in Disordered Granular Matter 
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The confining pressure P is perhaps the most important parameter controlling the properties 
of granular matter. Strongly compressed granular media are, in many respects, simple solids in 
which elastic perturbations travel as ordinary phonons. However, the speed of sound in granular 
aggregates continuously decreases as the confining pressure decreases, completely vanishing at the 
jamming-unjamming transition. This anomalous behavior suggests that the transport of energy at 
low pressures should not be dominated by phonons. In this work we use simulations and theory to 
show how the response of granular systems becomes increasingly nonlinear as pressure decreases. 
In the low pressure regime the elastic energy is found to be mainly transported through nonlinear 
waves and shocks. We numerically characterize the propagation speed, shape, and stability of these 
shocks, and model the dependence of the shock speed on pressure and impact intensity by a simple 
analytical approach. 



I. INTRODUCTION 

Granular systems are materials formed by conglomer- 
ates of macroscopic particles repelling through their mu- 
tual contact These, a priori, simple systems exhibit 
complex collective behavior and physical response. For 
example, depending on external conditions, granular me- 
dia are found in fluid or solid phases, with defined char- 
acteristic features [||. Despite the visual analogy with 
more conventional matter, the nature of these granular 
phases is far from being well understood. 

Key to better understanding the behavior and response 
of granular systems is the study of the mechanisms of en- 
ergy transport and dissipation. Experiments on clastic 
wave propagation have shown that the confining pres- 
sure P controls the speed of sound c of the medium 0]- 
Q, causing it to scale all the way to zero at low pressure. 
Several works have shown that for high confining pres- 
sures the speed of sound scales as c - P 1/e , a relation 
that can be rationalized from an effective medium the- 
ory. However, at low confining pressures a much faster 
power law is observed c ~ P 1 / 4 , whose physical origin 
still remains unclear. In general, the complexity of wave 
propagation in granular matter is evidenced by the coda- 
like signal spectrum of elastic waves, originated from the 
scattering of sound with the heterogeneities of the media 

From the theoretical point of view, much of the cur- 
rent understanding of granular mechanics comes from a 
simple microscopic approach, where granular media are 
modeled as amorphous packings of soft repulsive spheres 
0-@ ■ Using this simple model, the elastic response has 
been mainly studied in the linear regime, by applying 
a conventional normal mode expansion of the potential 
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energy [l(|-[2Tj]. Following this normal mode approxima- 
tion, different studies revealed that strongly compressed 
granular media have a vibrational spectrum which fol- 
lows the Debye law, and therefore, their elastic behavior 
parallels the one of simple atomic solids. 

Such studies also revealed an anomalous behav- 
ior in the elastic response as the pressure on the 
medium decreases, and the system approaches the jam- 
ming/unjamming transition. In this low pressure regime 
it was found that the vibrational modes resemble ordi- 
nary phonons only below a characteristic frequency w*, 
that vanishes as P goes to zero fl2l . Il3j |. Above uj* the 
modes are still extended but strongly scattered by disor- 
der, hence their ability to transport energy is impaired 
fl7l [20j . This suggests that granular media behave like 
a continuum elastic medium only above a characteristic 
length scale I* (thought to be related with u*) which 
diverges as P goes to zero [l3j . 

However, there is still another source of anomalous be- 
havior in the low pressure regime, which is universally 
displayed by granular media, but it has received far less 
attention. As a direct consequence of the nonlinear de- 
pendence of the local contact force on the grain defor- 
mations, the sound speed of a granular medium vanishes 
as P goes to zero (22|, [23|. This fact suggests that the 
transport of energy in the low pressure regime should not 
be dominated by phonons. Actually, linear sound does 
not even exist at the jamming point, and a phononic de- 
scription here looses its validity. Thus, we can expect 
that at low pressures the transport of energy in granu- 
lar media would be dominated by supersonic non-linear 
waves and shocks. 

This peculiar property of a vanishing sound speed 
at non-vanishing densities, not found in other systems 
like gases, liquids or solids, was first studied by V. 
Nesterenko, who conceived the name sonic vacuum for 
a media whose sound speed is exactly zero (22[, [HI- 
Nesterenko and coworkers mainly analyzed the conse- 
quences of the sonic vacuum in the transport of energy 
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in ordered granular chains [22j-[28j]. In such cases clastic 
energy was found to be transported in the form of non- 
linear solitary waves (characteristic energy pulses that 
propagate without shape distortion). 

Despite the results for ordered chains, the mechanisms 
of non-linear energy transport in disordered granular 
matter is largely unexplored [24|, Hf|- In particular, it is 
not clear what role the structural disorder of the medium 
plays. Moreover, no systematic study exist of the ex- 
pected crossover between the linear waves obtained at 
strong confining pressures and the nonlinear waves ex- 
pected to dominate in the low pressure regime, when 
the granular medium is close to the jamming/ unjamming 
transition. 

In this work, we study the non-linear transport of en- 
ergy in 2D disordered granular matter by using numer- 
ical simulations and theoretical calculations. Here we 
focus on the full response of two-dimensional granular 
packings when they are dynamically compressed by an 
external piston. These experiments cause the formation 
of steadily propagating fronts, without attenuation, nor 
deceleration, allowing the clear determination of the na- 
ture of the associated elastic waves (which can be linear, 
weakly non-linear, or strongly non- linear). By using this 
approach, we show how the elementary excitations of sys- 
tems at low pressures are shocks waves rather than ordi- 
nary linear phonons [29] . We provide a full characteriza- 
tion of these shocks, studying features like shape, prop- 
agation speed, and stability. We also present a simple 
model which accurately describes the propagation of the 
shocks, and the crossover from shocks to linear waves, as 
the pressure increases, or the driving strength decreases. 

The paper is organized as follow: In section II, we de- 
scribe the numerical model and the molecular dynamics 
simulations. In section III, we present the numerical re- 
sults related to the formation of shock waves by piston 
compression experiments. In section IV, we develop an 
analytical model which captures the speed of the shocks 
from the nearly linear to the strongly nonlinear regime 
using conservation laws. The appendices present an al- 
ternative approach to this behavior: a continuum equa- 
tion of motion (a simplification of Nesterenko's) is de- 
rived in Appendix A while its soliton and shock solutions 
are presented in Appendix B. 



II. NUMERICAL MODEL AND SIMULATIONS 



Materials with granularity, as diverse as sand, rocks, 
ceramics, foams, or emulsion, can exhibit a universal elas- 
tic response which is well described by a simple micro- 
scopic approach In all these systems the elasticity 
is governed by the soft contact interaction between the 
grains. In the simplest model of frictionlcss spherical 
grains, the contact potential between two grains i and j, 



located at positions Xi and xj , is written as jl5| : 

^• = ^% Q if %>0 (1) 
Uij = if 5ij < { ' 

where the overlap between particles is given by Sij = 
Ri + Rj — \xi — Xj\, where Ri is the radius of the i th 
particle. The exponent a gives the strength of the inter- 
action. The case a = 5/2 corresponds to the widely used 
Hertz's law. The interaction parameter = f r^+r ^tj 
is expressed in terms of the effective Young's modulus of 
the two particles, EL; see Ref. [l5[ for more details. 

We model granular media by means of jammed amor- 
phous packings of particles interacting through the 
Hertzian law (a = 5/2). The packings are prepared at 
a fixed pressure P, or equivalently, an average particle 
overlap <5jn~ P 2 / 3 , by using a molecular dynamics al- 
gorithm [15(. Similar structures are also commonly ob- 
tained through conjugate gradient algorithms. In order 
to avoid crystallization, the particle radii are uniformly 
distributed between 0.8 and 1.2 times their average ra- 
dius. Typical packings arc in the range of 10 3 to 10 4 
grains with various width to length ratios. 

In this work, we study the transport of energy by 
means of piston compression experiments. Here an ini- 
tially jammed configuration is continuously compressed 
by a piston which moves with a constant velocity up 
in the X direction, see the schematic in Fig. [1] The 
subsequent motion of the particles is obtained by numer- 
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FIG. 1. Piston compression experiments. A continuous com- 
pression of the system by a piston moving at velocity up leads 
to the formation of a pressure front (shock wave) travelling 
at a velocity vs- The shock compression process can be seen 
as a fast transition between two states (initial and shocked) 
of the granular system. 
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ical integration of Newton's equations of motion using 
the Verlet algorithm [30|]. We employ periodic boundary 
conditions in the Y direction and a fixed boundary on 
the right edge of the system. In our simulations, the unit 
of mass is set by fixing the grain density to unity. The 
effective particle Young modulus E* is set to one, which 
becomes the pressure unit. These choices ensure that the 
speed of sound inside the grain, v g is one (l5j . Lengths 
are measured in units of average particle diameter. 

Piston compression experiments have been previously 
used in studies of shock waves travelling through differ- 
ent media, like gases, liquids, or crystals [3l|-[33. In all 
these systems a travelling pressure front is observed as a 
consequence of the motion of the piston. In general, the 
velocity and shape of the front reveals important infor- 
mation about the mechanisms of energy transport. 



III. NUMERICAL RESULTS 
A. Phenomenology of Shocks 

Simulations at different compression velocities show 
that piston compression experiments induce the forma- 
tion of travelling pressure fronts (Fig. Q]). Ahead of 
the front there is an undistortcd static region where the 
particles remain at rest having their initial overlap Sq. 
This initial state is characterized by the original values 
of pressure Po , density po , strain 70 , and speed of sound 
Co- Behind the front we find a compressed region where 
the particles have acquired some kinetic energy. 

During their propagation the fronts remain essentially 
fiat, such that to a good approximation the compression 
experiments can be described as one dimensional pro- 
cesses. We present a detailed analysis and discussion of 
front stability and roughness at the end of the numerical 
results. 

Taking advantage of the flatness of the front, we av- 
erage over the transverse direction y, and obtain profiles 
of the averaged pressure P(x), density p(x), and the x 
component of the particle velocity u(x), as a function of 
the coordinate x. Figures 2a and 2b show the advance of 
a typical front, as seen through averaged profiles of pres- 
sure and velocity. Note that the velocity profiles show 
that behind the shock front the particles move on aver- 
age at the piston velocity up. 



In the steady state regime, the compression front at- 
tains a stationary shape described by an empirical fit, 
drawn as continuous lines in Fig. [2^, and[2jD. In general, 
the speed of propagation and the shape of the profiles de- 
pend on (1) the degree of non-linearity of the shocks and 
(2) the compression velocity up. In the next sections, 
we present a full characterization of such dependencies. 
By tracking the position of the profiles we find that af- 
ter transients have died out, the front travels linearly in 
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FIG. 2. Phenomelology of shock waves. Propagation of fronts 
as seen through the profiles of the averaged pressure P (a), 
and averaged particle velocity u normalized with the piston 
velocity up (b). The front moves in the direction of the 
arrows. Different curves corresponds to profiles at different 
times and the red lines are phenomenological fits, c) Position 
of the front as a function of time. Note that at long times the 
front moves linearly with time (red line), with a propagation 
speed vs- 



time t, propagating thus with a characteristic constant 
speed vs, see Fig. 3. The determination of the propaga- 
tion velocity allows the identification of the compression 
fronts with non-linear supersonic shock waves (vs cq), 
or linear waves (vs = Co). 

From the averaged profiles it is clear that the particles 
behind the front remain in a rather homogeneous state 
(see Fig. [3Ji and [2b)- This region is well characterized by 
the average values Ps, ps, 7s and cs- Thus, the passage 
of the front can be seen as a transition in the granular 
media from an initial (Pq, pq, 70, cq) to a compressed 
state (Ps, ps, 7s, cs). This behavior is strikingly dif- 
ferent from the profiles observed in crystals or crystalline 
arrays of grains, where large oscillations in pressure, den- 
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FIG. 3. Shock waves in a granular crystal, a) Cartoons 
of the pressure oscillations as a consequence of the in-phase 
motion of crystalline planes (red colors indicate regions of 
higher pressure), b) Propagation of the shock wave as seen 
through the particle velocity profiles. The front travels from 
left to right, c) Shock profiles in the presence of viscosity, for 
viscosity coefficients above (thick black line) and below (thin 
grey line) critical value. 



sity, or particle velocity are observed behind the fronts 

As an example, Fig. 3a shows the propagation of a 
shock wave obtained though the uniform compression of 
a granular crystal of hexagonal symmetry with zero pre- 
comprcssion (So = 0). Figure [3J5 shows the propagation 
of the shock as seen through the particle velocity profiles. 
Here the large coherent oscillations are caused by the in- 
phase motion of the crystalline planes. Note also that 
these oscillations change with time (compare grey and 
black curves corresponding to the time evolution), such 
that the profile never reaches a steady state. 

We believe that the underlying mechanism of equilibra- 
tion producing the homogenization of the profiles in the 
amorphous systems is intrinsically related with the struc- 
tural disorder. In our amorphous structures the disorder 
prevents any coherent oscillation, such that the peaks 
are washed out, leaving an homogeneous state. Similar 
homogenous states are also commonly observed in the 
compression of other unstructured systems, like gases or 
liquids [Hj]-[33|]. 

In ordered systems, homogeneous shock profiles are 
obtained when some source of dissipation is present 
[Hi , [13 ■ Figure shows a steady shock wave profile 
(black line) obtained when a viscous dissipation term, 



based on the relative velocities between particles, is 
added to the force: fij = rj\vi — Vj\, where r\ is the vis- 
cous coefficient. In this figure the viscosity coefficient has 
been chosen high enough to obtain a steadily propagat- 
ing profile without oscillations. For values of the viscos- 
ity coefficient smaller than a critical value, the obtained 
fronts reach steady profiles that still display some resid- 
ual oscillations (Fig. gray line). These conclusions, 
concerning the effect of viscosity on the shock profile, are 
demonstrated mathematically in Appendix B. 

Although our disordered granular structures do not 
have any source of real dissipation, the obtained uniform 
shock profiles suggest that the disorder of the system 
could be playing the role of an effective viscosity [29[ . 

B. Jump Conditions 

Mathematically, the propagation of fronts can be stud- 
ied through the conservation of mass, momentum, and 
energy: 

Pt + (pu) x = 
< (pu) t + (pv? + P) x = (2) 
^ (\pu 2 + U) t + [u(\pu 2 + U + P)] x = 

In the case of the propagation of uniform shocks, as found 
in our simulations (Fig. 2), the (approximate) uniformity 
of p, P and u ahead and behind the shock front allows 
the immediate integration of the equations, leading to 
a set of relationships known as the Rankine-Hugoniot 
jump conditions (3~l|-[33j: 

Povs = Ps{vs - up) 
^ po v 2 s + P = p s (up - v s ) 2 + P s 
ipo^l + Vs U + v s Pq = 

^Ps(vs - up) 3 + (v s - up) U s + (vs - up) P s 

(3) 

Here the results are expressed in the reference frame of 
the shock front. In this frame, the shock is in a steady 
state and the net flow of mass, momentum, and en- 
ergy toward the shock has to be zero. The quantities 
po, ps, Pq, Ps, and Uo,Us are respectively the density, 
pressure, and internal energy density before and after 
the passage of the shock. 

In the case of a Id chain, the density takes the simple 
form p = 2R-5 ana - t ne conservation of mass leads to 
a relation between the characteristic velocities up and 
vs, through the average radius of the particles, R, and 
the compressions in and ahead of the shock, 5s and So 
respectively: 

Vs = up — . (4) 

os - o 

Note that since the particle compression <5g is typically 
much smaller than its diameter 2R, Eq. Q implies that 
vs 3> up. In the limiting case of hard incompressible 
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spheres v$ — > oo. In reality, there is an upper cut-off to 
the speed v s given by the sound speed within the grain, 
that is equal in the case of steel to approximately 6000 
m/s . 

C. Shock's speed 

Commonly, the nature of the traveling fronts (linear or 
non-linear) depends on both piston velocity up and ini- 
tial pressure applied on the system Pq. Here, we analyze 
changes in the features of wave propagation by tracking 
variations in the front speed as a function of piston veloc- 
ity and the initial confining pressure Pq, which measures 
the distance of the uncompressed system to the unjam- 
ming transition, which occurs when Pq — > 0. This allows 
us to determine when nonlinear sound propagation takes 
over from the linear propagation. 

Figure [4ji shows the results of this study. Note that 
for any finite initial pressure Po, at low compression ve- 
locities up the fronts travel with velocities vs which are 
roughly independent of up . This is evidence of a regime 
of linear/ weakly nonlinear wave propagation. Here the 
fronts travel at approximately the longitudinal speed of 
sound of the media v$ ~ cq. In this regime, the propa- 
gation speed is simply controlled by the initial confining 
pressure Po or initial overlap Sq. 

The inset of Fig. [4ji shows the obtained variation of 
the speed of sound cq with the initial overlap Sq- The 
functional variation of cq with Sq can be rationalized us- 
ing scaling arguments. Note that in general cq ~ \f~B~, 
where the bulk modulus B = and P = |f . The 
change in volume dV scales linearly with Sq, the aver- 
age overlap between particles, while the energy scales as 
E ~ Sq, see Eq. [TJ Upon setting a = 5/2, we obtain the 
pressure dependence of the longitudinal speed of sound 
Co ~ S\/ A ~ P 1 / 6 valid for Hertzian interactions [l5j |. 
The inset of Fig. shows that the numerical data for 
Co, represented by symbols, is consistent with the Sq^ 4 
scaling, which is shown as a continuous line. 

As the compression velocity up is increased, the speed 
of the front gradually grows due to a departure from lin- 
earity in the mechanisms of energy transport (Fig. HJt). 
In this regime elastic energy starts to be mainly propa- 
gated through supersonic shock waves (vs > Co). 

Note that Fig. @Ji clearly shows that the different data 
for vs vs up has the same qualitative structure. This 
allows the creation of a single master curve, with the dif- 
ferent data collapsed, as shown in Fig. [4jx We achieve 
this upon the rescaling of the vs and up axis. The speed 
of propagation is normalized with the speed of sound of 
the medium M = vs/cq, where M is the Mach number. 
The up axis is rescaled by a pressure-dependent veloc- 
ity scale up, which marks the crossover between linear 
acoustic waves and shocks. 

Above up , the transport of energy starts to enter into 
a new regime dominated by strong shock waves, which 
travel at two or three times the speed of sound (Fig. 



Hp). In this regime we find a power law relation between 

particle and front velocities vs ~ u p /5 (or M ~ uj 5 ) as 
indicated by the dashed line. 

The main features of this fast compression regime 
(up up) can also be rationalized through scaling. 
Since up, R and Sq are all known, we need one addi- 
tional relation which combined with Eq. Q will make a 
definite prediction for the shock speed. We note that for 
strong shocks, the propagating front generates a charac- 
teristic compression S Sq and a corresponding increase 
in the kinetic energy. In general for strong shock waves 
the kinetic and potential e nerg ies in the compressed re- 
gion are of the same order (|32j). In our case this implies 
up ~ S 5 / 2 . We have tested numerically that this non- 
trivial proportionality relation exists for strong deforma- 
tions [291 ] . Upon combining the balance between kinetic 
and potential energy with Eq. (U), one readily obtains 

1 /5 

the observed power law vs ~ u p . 

We conclude that by controlling Sq or P, which param- 
eterize the distance to the jamming point (at P = and 
Sq =0), we can tune up and the onset of the strongly 
non-linear response of the packings. Our key numerical 
findings on the shock velocity summarized in Fig. U can 
be grasped from scaling near the jamming point. 




FIG. 4. Propagation speed, a) Speed of the front vs versus 
particle velocity up measured in units of v g , the sound speed 
within the grain, for decreasing particle overlap So- The in- 
set shows the variation of the sound speed Co with the initial 
overlap <5o (data) and the scaling Co ~ delta\ (line), b) Nor- 
malized plot obtained be rescaling vs with co vs(0) and up 
with Up. The dashed line indicates the power law vs ~ lip 1 / 5 
characteristic of a sonic vacuum. The inset shows the varia- 
tion of up with the initial overlap So (data) and the scaling 
up ~ 6^ 4 (line). 
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D. Shock's shape and width 

Figure [5] shows typical shock profiles observed for weak 
and strong shock waves. In general, there are two main 
significant changes in the profiles due to the increase in 
the non-linearity of the phenomena. First, the width of 
the shocks decreases when increasing the non-linearity. 
Second, while symmetric profiles are obtained for weak 
shocks, strong shocks are associated with much more 
asymmetric profiles. Here we found that in all these 
cases, the profile of the shocks can be well approximated 
through the phenomenological expression: 



u(X) 



1 + expfri {X - X )] 



(5) 



where u p is the compression velocity, Xq is the position 
of the shock wave, and the p^s are fitting parameters nec- 
essary to describe the profile [29} . While in general weak 
(symmetric) shock profiles are well described by only one 
parameter p\ , the description of strong shocks requires at 
least three fitting parameters Pi-p?,. The continuous lines 
of Fig. [5] shows typical fits for weak and strong shocks. 

In order to characterize the features of the shock pro- 
files as a function of non-linearity, we systematically mea- 
sure shock widths and symmetries for the different com- 
pression velocities, on systems with varying initial pres- 
sure. We define the shock's width A and asymmetry 
parameter Q by the expression (39j.(40j: 



(6) 



max(^) 



q = Q2^ u ( x ) dx 

® l ST 00 [up-u(n)]<te; 



where x* is such that u{x*) = u p /2. Figure |Hk shows 
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FIG. 5. Typical shock profiles for a) weak shock and b) 
strong shock. The continuous lines represent phenomenolog- 
ical fits. 
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FIG. 6. Shape characterization of the shock profiles, a) 
Scheme of the definition of shock's width and asymmetry. 
Variation of shock width (b) and asymmetry (c) with the 
Mach number for the different studied shocks. The color code 
is the same as in figure [4] 



a schematic of the shock profile illustrating the physical 
meaning of these definitions. Figures |6b and Eh show the 
results of shock's width and asymmetry, as a function of 
the Mach number of the shocks. The width and symme- 
try of the shocks monotonically decrease as the nonlin- 
earity increases. Note that while the width reaches 50 
particles for weak shocks, it reduces to approximately 10 
particles for strong shocks. Also it is clear that weak 
shocks are rather symmetric (Q ~ 1) and strong shocks 
are asymmetric Q < 1 due to the slow relaxation behind 
the front. 



E. Characterization of the shocked state 

We have seen how the passage of the pressure front 
produces a transition in the granular media from an ini- 
tial state to a different shocked state. Since the width of 
the shock waves is generally small, the shock compression 
process is usually considered to be a fast adiabatic tran- 
sition [Ilj]-[33j]; anows one to trace out the adiabatic 
equation of state by measuring the properties of shocks. 

Generally, the features of the compressed region are 
studied through the shock adiabat, or Hugoniot curve, 
which is the locus of final shocked states. The Hugoniot 
is an adiabatic line in the P—V plane (see Fig. [7^). Each 
point in this curve indicates the volume and pressure in 
the medium in the final shocked states. Note that during 
the compression experiment, which induces the transi- 
tion from an initial state (Vo,Po) before the shock to a 
final state (Vg, Ps) after the shock has passed, the system 
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a) b) 




FIG. 7. a) Schematic of a general shock Hugoniot curve 
(shock adiabat). b) Numerically obtained shock Hugoniot for 
systems with different initial pressures as obtained from the 
jump conditions. Dashed lines are guides to the eyes. 



does not follow the Hugoniot curve in the intermediate 
states. In general the evolution is given by the Rayleigh 
line which is the straight line connecting initial and final 
states in the P — V plane (see Fig. [T^t) . 

As the shock compression process is adiabatic, the 
Hugoniot curve is often used to construct the equation 
of state of a given material. Following this idea, shock 
compression experiments have been designed and used for 
years to study the equation of state of solids and minerals 
at high pressures |33l |. 

In our case the shock Hugoniot can be directly ob- 
tained from the simulations (by direct inspection of the 
final shocked states), or by using the Rankine- Hugoniot 
jump conditions. The conservation of mass and momen- 
tum across the shock front Eq. [3J can be rewritten re- 
lating the pressure and volume in the material with the 
particle and front velocities (3l|-[33|: 

V/V = 1 - up/vs 
P - P = p up v s 

When the front velocity is known as a function of particle 
velocity v$ = vs{up), these relations become an implicit 
expression for the Hugoniot curve. Figure [7}d shows the 
shock Hugoniot curve as obtained by using the Rankine- 
Hugoniot conditions. In what follows, we use this numer- 
ically determined Hugoniot adiabat to study the stability 
of the shock waves for the case of small disturbances in 
the shape of the shock front. 

F. Shock's stability, roughness, and focusing 

The stability of shock waves against shape perturba- 
tions has been one of the interesting and important topics 
of compressible flow. The theoretical investigations on 
this topic started with the seminal work of Dyakov [4l[ , 



based on a normal mode expansion of the front pertur- 
bations, later re-examined and presented in more detail 
by Swan and Fowles [42| . 

Dyakov performed a linear stability analysis obtaining 
conditions for shock stability under which a shock wave 
remains stable, and front perturbations decreases expo- 
nentially in time. In this analysis the stability of the 
shocks is related with the slope of the Hugoniot adiabat, 
such that a shock wave remains stable whenever these 
conditions hold HJ, 

3 2 dV/dP < 1 + 2M S 
3 2 dV/dP > -1 1 ' 

In these expressions the derivatives are taken along the 
Hugoniot adiabat, Ms is the Mach number of the shock 
respect to the material behind the front, and j 2 is the 
slope of the Rayleigh line. In general, when the lower 
limit of the above equation is not satisfied the shock 
wave splits into two waves travelling in the same direc- 
tion (splitting instability). The upper limit corresponds 
to the corrugation instability, where front modulations 
are exponentially amplified. 

Note the simplicity and generality of this result. Once 
the Hugoniot curve is known these conditions are gen- 
eral and can be applied to any material. By using our 
numerically obtained shock Hugoniot, we can perform 
the Dyakov test stability analysis. Figure [8] shows our 
results, which suggest that the shock waves studied in 
this work are stable for the broad range of compressions 
and impact intensity analyzed. 

Although the shock waves are stable, the disorder of 
the granular packings induce the roughening of the front 
(compare the fronts of Fig. Q]with Fig. [3Ji). In general, 
in crystalline solids shock front's may also become irreg- 
ular when travelling through poly-crystalline structures. 
Although not a proper instability, changes in the speed of 
sound with the lattice orientation can induce the corru- 
gation and roughening of the shock front [43| , [H| . In the 
amorphous structures considered here, the complete lack 




FIG. 8. Dyakov's stability analysis. Note that all the shock 
waves belong to the stable region of the phase diagram. 
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time 

FIG. 9. Evolution of the front's roughess as a function of time 
for different Mach numbers. Note that when the nonlinearity 
of the wave is large, the roughness reaches a plateau at a small 
value. 



of symmetry produces a rather isotropic speed of sound 
such that irregularities in the front are restricted to a few 
particles diameters. 

In order to analyze variation in front roughnesses we 
define the average front position (x) and roughness Ax 
as: 

(x) = {x (y)) , , 

A x =((x (y)-(x)) 2 ) 1/2 I > 

where xo(y) is the position of the front at y (sec Fig. 
Figure |9}d shows the time evolution of the front rough- 
ness for shocks of varying Mach number. Interestingly, 
although the roughness of the fronts seems to reach max- 
imum values less than 5 — 6 particles, the value of the 
plateaus clearly depends on the nonlinearity of the wave. 
In particular, our simulations suggest that the more pro- 
nounced the nonlinearity of the wave is, the less rough 
its front is. 



IV. ANALYTICAL THEORY OF THE SHOCK 
SPEED 

In this section, we analytically describe the relation 
between the shock speed and strength, shown in Fig. |H 
using the jump conditions Eq. [3l Appendix B gives 
an alternative derivation. The first two equations can 
be used to predict the speed of the shock if we assume 
that the pressure is given entirely by the elastic forces 
between the particles. We describe the system with a 
one-dimensional one; then any interface cuts only one 
pair of particles from each other, and the pressure is just 
identical to the force between these two particles, and 
is given by Hertz's law, P = e<5 a_1 . Combining this 
equation of state with the expression for the density p = 
2ft_$ 7 the mass and momentum conservation conditions 
(first two relations in Eq. [3]) imply: 

c y a - 1 Ss/5o - 1 

where Co represents the speed of sound before the passage 
of the shock wave. Together, Eqs. @] and [TT] can be 
seen as a parametric relation between front and particle 
velocities, where the overlap 5s produced by the passage 
of the front is the parameter. Such a parametric plot 
of vs versus up is drawn as a continuous line on the 
numerical data in Fig. |4Jd. This comparison shows that 
Eqs. 2] and [TTJ are in excellent agreement with the results 
of our numerical experiments on shock propagation. 

Despite the good agreement with the numerical results, 
the third jump condition is not satisfied if we assume the 
previous expressions for Vs, up, and energy density per 
unit length U = (5 a -j^ according to Hertz's law. In the 
simulations the particles also have a random motion su- 
perimposed on the average velocity up. This produces 
an excess in the pressure and energy of the particles. In 
this case, the internal energy has an additional contri- 
bution, which can restore conservation of energy. More 
precisely, the random motion is described by an effective 
temperature Tg, and hence the speed of the shock is no 
longer overdetermined by the three jump conditions be- 
cause the total number of unknowns, r>s, 5 s and Tg is 
now equal to the number of conditions. 

These heating effects should change the shock speed, 
but not drastically. The prediction above was based en- 
tirely on conservation of mass and momentum, which 
must still hold. However, the pressure will have a contri- 
bution from the thermal motion, and this can affect the 
conclusion. To estimate how big of an effect the thermal 
motion has, let us estimate how much energy is dissi- 
pated. First, the velocities may be eliminated from the 
three conservation laws to give: 

PqUs ~ PsUq = Pq + Ps ^ 
ps - po 2 

(this is obtained by solving for the velocities from the 
first two conservation laws and substituting into the 
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third). Taking Pq, P$ and Uq to have approximately 
the Hertzian form (i.e., neglecting the temperature), we 
can solve for Us- The excess over the Hertzian potential 
energy, p ffjf is then the thermal energy, which can be 
expressed entirely in terms of So and 5s ■ 

U thermal — r\ r-» c 

2R-6 S 

(13) 

When a is close to 2, this energy excess becomes insignif- 
icant since this expression approaches 0. In particular, 
when So = and a = §, L thermal is — § of the po- 
tential energy, explaining why our approximation agrees 
well with the numerics. 

To summarize, the speed of a shock as a function of 
the compression can be estimated from conservation laws 
if one neglects the heating of the particles. In the limit 
where a is close to 2 (as in the Hertzian case) this ap- 
proximation is quite accurate. 

One could go beyond this approximation and deter- 
mine the exact relation between speed and compression 
in a shock. First the thermal contribution to both pres- 
sure and energy has to be predicted as a function of tem- 
perature. All the parameters of the shock can then be 
determined by solving all three jump conditions simulta- 
neously. 



V. CONCLUSIONS 

In this work, we have used simulations of piston com- 
pression experiments to unveil the mechanism of energy 
transport in granular systems composed of soft friction- 
less spheres. While at high confining pressures the elas- 
tic response of the granular system is basically linear, 
the low pressure regime is found to be dominated by 
non-linearities. Thus, for low confining pressures, when 
the system is close to the jamming-unjamming transition, 
the basic elastic excitations are supersonic shock waves, 
rather than linear phonons. 

We have presented a detailed characterization of the 
weakly and strongly non-linear shocks, studying their 
propagation speed, width, shape, and stability. It is in- 
teresting to note that the propagation speed of shocks 
does not depend on friction or dissipation. Thus, experi- 
ments should focus on this quantity which is only affected 
by the contact interaction law. In general, the presence 
of friction or dissipation will affect other features like the 
width of shocks. 

Finally, we developed a mathematical model that cor- 
rectly captures the universal nonlinear response of granu- 
lar packings around the jamming-unjamming transition. 
In particular, we have used this model to describe the 
dependence of the propagation speed on pressure and 
driving. 

We conclude by saying that the elastic response of 
granular matter present unique features, not found in 
other condensed matter systems like solids, liquids, or 



gases. The athcrmal character of these systems allows the 
attainment of a state with zero pressure at non- vanishing 
density, the jamming point, which corresponds to a sonic 
vacuum. In this state the response of the system is char- 
acterized by a complete lack of linear behavior. Here we 
have shown how this fully nonlinear state also controls 
the elastic response of granular system in the low pres- 
sure regime. 
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Appendix A: Equation of motion 

The Lagrangian of a Id chain composed of soft spheres 
with repulsive interactions, as in Eq. 1, reads: 

£ = E^*n-%n-^ + l)", (Al) 

* — ' 2 a 

n 

where ip n is the displacement of the n particle, m is the 
mass, e the interaction's constant, and a the interaction's 
exponent (a — 5/2 for Hertzian interactions). The con- 
tinuum limit of the system is obtained by considering the 
limit where the radius of the spheres R — > 0. This limit 
is justified when the wavelength of waves is large com- 
pared to the size of the spheres. (As Nesterenko found, 
the wavelength is actually a fixed multiple of the size of 
the spheres, but the multiple is large enough that the 
continuum approximation is justified.) Here we proceed 
to replace the ip n 's by a continuous function and perform 
the continuum approximation in the Lagrangian, before 
deriving the equation of motion. 

We first define a continuous function <p interpolating 
through the ij) n '8. Although the most obvious choice is 
4>(2Rn) = ip n , it is not the easiest to implement. In 
this case we have i/) n — ip n +i = —2R(f>'{x) — 2R 2 <f)"(x) — 
. . . where x = 2Rn is the equilibrium position of the 
n th sphere. Thus, this approximation will complicate the 
calculation of the potential energy, because we need to 
use the binomial theorem to raise tp n — ip n +i to the a th 
power. 

In order to get a simpler equation we take the contin- 
uum limit differently Define (f> by prescribing its deriva- 
tive instead: 

V>„+i -il) n = 2R<t>'{2R{n + 1)) (A2) 

In this case the calculation of the potential energy will be 
simple, while the kinetic energy will receive corrections 
from an expansion in powers of R. (In Nesterenko's ver- 
sion of the equation, dispersion effects appeared in the 
potential term so they were nonlinear and more compli- 
cated.) We first invert Eq. (|A2[) in order to express ?/> n 
in terms of <j>. To do that we make a Taylor expansion of 
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the left-hand- side about n + |. This gives (up to terms 
of order (2R) 4 ): 

1/(2R{n+±)) + MV(2ii(n + \)) = <j>>(2R(n + ±)). 

(A3) 

Integration of both sides with respect to x leads to 4>(x) = 
(1 + Solving for ijj in terms of 4> gives: 



ip(x) 



1 



R 2 d 2 
6 dx' 



■<f>(x). 



(A4) 



Expanding the ratio in a geometric series leads to ip( x ) ~ 
4>{x) — ^?-(f)"(x) .... Thus the kinetic energy term can be 



approximated by: 
1 



2 - ^(*) 2 



<j)(ri)<j)"( 



(A5) 



The approximated continuum Lagrangian is obtained 
by replacing the tf> in terms of <fi and the sum by an 
integral: 





h" E(2R) °( *t 

ma , 




T 6 


' 2# 




2 £ (2i?)« \ 

771 CK / 


dx 




2R' 



(A6) 

Finally, by using the Euler-Lagrange equation, we obtain 
the equation of motion: 



3 







(A7) 



In terms of the compression field 5 = —2R<f>' this equation 
takes the form: 



d-—S"-^^[8 a - 1 }" = (A8) 
3 m 

This equation is simpler to the equation originally de- 
rived by Nesterenko; there is more than one way to take 
the continuum limit when the wavelength is the same 
order as the particle spacing. Here we found that this 
simpler equation is enough to describe the simulation re- 
sults. 




FIG. 10. Effective potential surface. For vs > c the ini- 
tial state becomes a local maxima allowing the description of 
travelling perturbations. 



Integrating twice and applying the boundary conditions 
S(i —> — oo) = 5 , S(t —> — oo) = 5(i —> — oo) = 0, leads 
to the following expression: 



R 2 ■■ 4eR 2 
6v s 



AeR 2 

T d Q 

mvi 



This last equation can be rewritten in the form: 



#8 _ 
dt 2 ~ 

where the field W is given by 
12e 



dW 
~dT' 



(Bl) 



3-1 



W{5,v s ) = — (<$*- 5«) ^(S-Sof 
ma ZR Z 

-— 6^(6 -6 ). 



(B2) 



Thus, the problem of traveling steady waves is mapped to 
the motion of an effective particle moving (in time t) in 
a potential field W. Note that in this case the potential 
not only depends on the compression 5, but also on the 
propagation speed vs- 



Appendix B: Steady state solutions 

The travelling of stationary perturbations having a 
propagation velocity vs can be considered by replacing 
the ansatz S(x,t) = 5(t = t — ^L) in the equation of 

motion. This function 5(t) describes what an observer 
at a fixed point of space sees as the shock passes by. It 
satisfies the ordinary differential equation: 

R 2 d 4 S 9 dS e d 2 ^"' 1 } 
3tj| dt A s dt mv% dt 2 



We now use this analogy to study the properties of 
travelling waves. Figure [10] shows the shape of the po- 
tential as a function of compression and propagation 
speed. In the initial state, the system is undistorted 
(S = So) and the effective particle is located at the ex- 
treme W(5o,v s ) = 0, at rest. For the cases where this 
initial state is a local minima, the effective particle re- 
mains at this initial position for all times. This means 
that in these cases steady perturbations cannot travel 
through the system. Thus, in order to describe travel- 
ling perturbations the initial state needs to be a local 
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maximum such that Wss{6q,vs) < 0, and the speed of 
propagation needs to satisfy: 



vi > 



4i? 2 (a- l)e, 



(B3) 



Physically, this condition means that travelling perturba- 
tions will propagate at the speed of sound of the medium 

a-2 

c = 2RyJ (a — l)e/mS 2 , or faster. Figure [TU] shows 
the change in the nature of the initial state, from a local 
minima to a local maxima, as the speed of propagation 
of the waves is increased. 



a) 




1. Solitary waves 

Now consider a travelling perturbation which satisfies 
the above condition (vs > cq). Figure ITTb shows the 
shape of the potential W(6) for a given front velocity vs- 
Far away, the undistorted system corresponds to the ef- 
fective particle located at the maxima of the potential 
W = 0. The distortion induced by the passage of the 
wave can be seen as the particle rolling down the poten- 
tial. The effective particle will move up to the turning 
point, and then the conservation of energy implies that 
the particle would roll back to the initial state. 

This motion of the effective particle corresponds to a 
solitary wave travelling trough the system with velocity 
vs (Fig. [TTb). Note that these solitary waves are the 
only steady state solutions of the system. The speed 
of propagation can be obtained from the turning point 
condition W(6sol,vs) = 0, where 5 sol is the maximum 
compression due to the passage of the wave: 



Co 



a(a — 1) 



u SOL °0 



aS 



S%- 2 (5sol-6 ) 2 S sol 



(B4) 



2. Shock waves 

The previous analysis shows that in the absence of any 
source of dissipation the motion of the effective particle 
conserves its energy, and the only steadily propagating 
solutions are solitary waves. In that case homogeneous 
shocks, as found through the simulations in this work, 
are not possible solutions of the equation of motion. 

In order to obtain steady state shock solutions, some 
sort of dissipation is needed. In general, the effect of dis- 
sipation can be taken into account by adding new terms 
in equation IA71 For example, if a dissipation term (due 
to relative motion of the spheres) is considered [37|, the 
equation of motion of the system will take the form: 



3 



4R 2 s 



[(-* 



Ma — li ' 



4i? 2 ?y0" = 0, (B5) 



b) 



solitary wave 

SL 



c) 



shock wave 





FIG. 11. Steady state solutions of the non-linear wave equa- 
tion, a) Potential well for a given propagation speed vs- b) 
Solitary wave solutions obtained in te absence of viscosity, c) 
Shock wave solutions obtained for viscous dissipation. 



where 77 is a damping coefficient. 

Now if we look for steady state travelling perturba- 
tions with the same procedure as before, we obtain the 
following equation: 



dl 2 



dW dS 

'Is" 127 ] £ 



(B6) 



Here the potential W is the same as obtained before 
(equation (|B2|I ). Thus, for this particular form of dis- 
sipation, we have explicitly shown that the problem of 
steadily propagating waves is mapped to the problem of 
an effective particle moving in a potential, but now the 
particle also moves under the effect of friction. 

In this case, the presence of friction will take out the 
energy of the effective particle, inducing the oscillation 
to decay to the minimum of the potential well. These 
solutions correspond to oscillatory shock waves (Fig. [TTb . 
black line). 

If the viscosity is high enough, the motion of the effec- 
tive particle will become overdamped, and the particle 
directly moves to the minimum of the potential, without 
performing any oscillation. This motion corresponds to 
traveling regular shock profiles (Fig. [TTb, green line). 

For the shock waves the relation between propagation 

velocity and induced compression can be obtained by the 

condition 4§- c = 0, which corresponds to the mini- 
ad I 8s,vs 

mum of the potential. This condition leads to: 



vs_ 

CO 



1 {Ss/SoY 



- 1 



a — 1 



Ss/So 



(B7) 
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where 5s is the maximum compression after the shock 
wave has passed. In Sec. IV, imposing two of the jump 
conditions (conservation of momentum and mass, but not 
energy) produced the same result. The dissipation term 
we have added preserves the same two conservation laws, 



explaining the agreement. In Sec. IV, we justified the 
disappearance of energy not by introducing friction (since 
the simulations did not include it) but by assuming that 
the energy is converted to random motion of the particles. 
A more accurate treatment of this hypothesis would have 
to account for some additional thermal pressure. 
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